%% Import Results from Fortran into Matlab
% The Gender Investment Gap over the Life-Cycle
% Annika Bacher, annika.bacher@bi.no
% This version: May 2024

%% housekeeping !!! SPECIFY PATHS HERE !!!

        clear all;
        clc;
        close all;

% path to exported results from fortran
cd('');

% path to data targets 
datafolder = '';

% path to store figures 
figurefolder = '';

%% specify array lengths

        JJ     = 36;        % working age periods
        JR     = 20;        % retirement periods
        ncc    = 100;       % asset grid
        naa    = 2;         % education levels
        nzz    = 9;         % toal productivity states
        nz     = 3;         % persistent productivity states
        minage = 30;        % minimum age
        nrisk  = 20;        % risky asset grid 
        SS     = 25000;     % number of households to simulate
        
        age    = (minage-1)+linspace(1,JJ,JJ);
        cgrid  = linspace(1,ncc,ncc);
        prodgrid  = linspace(1,nzz,nzz);

%% read in life-cycle profiles from model

    % Stock Market Participation Rate
    fid = fopen('base\cohort\SMPcoh_sf.txt');
    SMPcoh_sf  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);  

    fid = fopen('base\cohort\SMPcoh_sm.txt');
    SMPcoh_sm  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);  
    
    fid = fopen('base\cohort\SMPcoh_c.txt');
    SMPcoh_c  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);
 
    % Assets
    fid = fopen('base\cohort\cashcoh_sf.txt');
    cashcoh_sf  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);  

    fid = fopen('base\cohort\cashcoh_sm.txt');
    cashcoh_sm  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);  
    
    fid = fopen('base\cohort\cashcoh_c.txt');
    cashcoh_c  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    % Unconditional Equity Share
    fid = fopen('base\cohort\ucsharecoh_sf.txt');
    ucsharecoh_sf  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);  

    fid = fopen('base\cohort\ucsharecoh_sm.txt');
    ucsharecoh_sm  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);  
    
    fid = fopen('base\cohort\ucsharecoh_c.txt');
    ucsharecoh_c  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    % Conditional Equity Share
    fid = fopen('base\cohort\sharecoh_sf.txt');
    sharecoh_sf  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);  

    fid = fopen('base\cohort\sharecoh_sm.txt');
    sharecoh_sm  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);  
    
    fid = fopen('base\cohort\sharecoh_c.txt');
    sharecoh_c  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

%% read in life-cycle targets from data

% Stock Market Participation Rates
SMPcoh_s_dat = csvread([datafolder '/SMP_single.csv']);
SMPcoh_sf_dat = SMPcoh_s_dat(JJ+1:end,3)';
SMPcoh_sm_dat = SMPcoh_s_dat(1:JJ,3)';

SMPcoh_c_dat = csvread([datafolder '/SMP_couple.csv']);
ageav = SMPcoh_c_dat(:,2);
SMP_c_dat = SMPcoh_c_dat(1:2:end,3);

SMPcoh_s_ci_dat = csvread([datafolder '/SMP_single_conf.csv']);
SMPcoh_sf_cl_dat = SMPcoh_s_ci_dat(JJ+1:end,3)';
SMPcoh_sf_ch_dat = SMPcoh_s_ci_dat(JJ+1:end,4)';
SMPcoh_sm_cl_dat = SMPcoh_s_ci_dat(1:JJ,3)';
SMPcoh_sm_ch_dat = SMPcoh_s_ci_dat(1:JJ,4)';

SMPcoh_c_ci_dat = csvread([datafolder '/SMP_couple_conf.csv']);
SMPcoh_c_cl_dat = SMPcoh_c_ci_dat(1:2:end,3);
SMPcoh_c_ch_dat = SMPcoh_c_ci_dat(1:2:end,4);

% Net Worth
nw_s_dat = csvread([datafolder '/mnw_single.csv']);
nw_sf_dat = nw_s_dat(JJ+1:end,3)';
nw_sm_dat = nw_s_dat(1:JJ,3)';
 
nw_c_dat = csvread([datafolder '/mnw_couple.csv']);
nw_c_dat = nw_c_dat(1:2:end,3);

nw_s_ci_dat = csvread([datafolder '/mnw_single_conf.csv']);
nw_sf_cl_dat = nw_s_ci_dat(JJ+1:end,3)';
nw_sf_ch_dat = nw_s_ci_dat(JJ+1:end,4)';
nw_sm_cl_dat = nw_s_ci_dat(1:JJ,3)';
nw_sm_ch_dat = nw_s_ci_dat(1:JJ,4)';

nw_c_ci_dat = csvread([datafolder '/mnw_couple_conf.csv']);
nw_c_cl_dat = nw_c_ci_dat(1:2:end,3);
nw_c_ch_dat = nw_c_ci_dat(1:2:end,4);

% Unconditional Equity Share
ucsharecoh_s_dat = csvread([datafolder '/share_single.csv']);
ucsharecoh_sf_dat = ucsharecoh_s_dat(JJ+1:end,3)';
ucsharecoh_sm_dat = ucsharecoh_s_dat(1:JJ,3)';

ucsharecoh_c_dat = csvread([datafolder '/share_couple.csv']);
ucsharecoh_c_dat = ucsharecoh_c_dat(1:2:end,3);

ucsharecoh_s_ci_dat = csvread([datafolder '/share_single_conf.csv']);
ucsharecoh_sf_cl_dat = ucsharecoh_s_ci_dat(JJ+1:end,3)';
ucsharecoh_sf_ch_dat = ucsharecoh_s_ci_dat(JJ+1:end,4)';
ucsharecoh_sm_cl_dat = ucsharecoh_s_ci_dat(1:JJ,3)';
ucsharecoh_sm_ch_dat = ucsharecoh_s_ci_dat(1:JJ,4)';

ucsharecoh_c_ci_dat = csvread([datafolder '/share_couple_conf.csv']);
ucsharecoh_c_cl_dat = ucsharecoh_c_ci_dat(1:2:end,3);
ucsharecoh_c_ch_dat = ucsharecoh_c_ci_dat(1:2:end,4);

% Conditional Risky Share
sharecoh_s_dat = csvread([datafolder '/share_an_single.csv']);
sharecoh_sf_dat = sharecoh_s_dat(JJ+1:end,3)';
sharecoh_sm_dat = sharecoh_s_dat(1:JJ,3)';

sharecoh_c_dat = csvread([datafolder '/share_an_couple.csv']);
sharecoh_c_dat = sharecoh_c_dat(1:2:end,3);

sharecoh_s_ci_dat = csvread([datafolder '/share_an_single_conf.csv']);
sharecoh_sf_cl_dat = sharecoh_s_ci_dat(JJ+1:end,3)';
sharecoh_sf_ch_dat = sharecoh_s_ci_dat(JJ+1:end,4)';
sharecoh_sm_cl_dat = sharecoh_s_ci_dat(1:JJ,3)';
sharecoh_sm_ch_dat = sharecoh_s_ci_dat(1:JJ,4)';

sharecoh_c_ci_dat = csvread([datafolder '/share_an_couple_conf.csv']);
sharecoh_c_cl_dat = sharecoh_c_ci_dat(1:2:end,3);
sharecoh_c_ch_dat = sharecoh_c_ci_dat(1:2:end,4);

 %%%% smooth net worth plot for singles %%%%
      %%%% moving average of 3 %%%%

   nw_sm_dat_smooth = movmean(nw_sm_dat,3);
   nw_sm_cl_sm = movmean(nw_sm_cl_dat,3);
   nw_sm_ch_sm = movmean(nw_sm_ch_dat,3);

   nw_sf_dat_smooth = movmean(nw_sf_dat,3);
   nw_sf_cl_sm = movmean(nw_sf_cl_dat,3);
   nw_sf_ch_sm = movmean(nw_sf_ch_dat,3);


%% for HC return: read in chains of shocks (productivity, marital, risky asset, aggregate state)

% NOTE: simulation is by gender (i.e., simulated for 25000 men and women each who
% may switch marital status)

 % risky asset  return (men and women)
 fid = fopen('hc_return\simdata\chain_retf.txt', 'r');
        chain_retf  = cell2mat(textscan(fid, '%f')).';
        fclose(fid);  
        chain_retf = reshape(chain_retf,SS,JJ+1);

 fid = fopen('hc_return\simdata\chain_retm.txt', 'r');
        chain_retm  = cell2mat(textscan(fid, '%f')).';
        fclose(fid);  
        chain_retm = reshape(chain_retm,SS,JJ+1);

 % marital status (men and women)  
 fid = fopen('hc_return\simdata\chain_martf.txt', 'r');
        chain_martf  = cell2mat(textscan(fid, '%f')).';
        fclose(fid);  
        chain_martf = reshape(chain_martf,SS,JJ+1);
 
 fid = fopen('hc_return\simdata\chain_martm.txt', 'r');
        chain_martm  = cell2mat(textscan(fid, '%f')).';
        fclose(fid);  
        chain_martm = reshape(chain_martm,SS,JJ+1);      
  
 % labor productivity (men, women, couples)
 fid = fopen('hc_return\simdata\chain_prodf.txt', 'r');
        chain_prodf  = cell2mat(textscan(fid, '%f')).';
        fclose(fid);  
        chain_prodf = reshape(chain_prodf,SS,JJ+1);
        
 fid = fopen('hc_return\simdata\chain_prodm.txt', 'r');
        chain_prodm  = cell2mat(textscan(fid, '%f')).';
        fclose(fid);  
        chain_prodm = reshape(chain_prodm,SS,JJ+1);

  fid = fopen('hc_return\simdata\chain_prodc.txt', 'r');
        chain_prodc  = cell2mat(textscan(fid, '%f')).';
        fclose(fid);  
        chain_prodc = reshape(chain_prodc,2*SS,JJ+1);

  % aggregate state
  fid = fopen('hc_return\simdata\chain_agg.txt', 'r');
        chain_agg  = cell2mat(textscan(fid, '%f')).';
        fclose(fid);  
        chain_agg = reshape(chain_agg,SS,JJ+1);
        
%% compute correlation between stock and HC return (Appendix C.1)
 
% NOTE: only done for singles!

   % import simulated human capital values (single women and single men)
   fid = fopen('hc_return\simdata\hcval_path_sf.txt', 'r');
        hcval_path_sf  = cell2mat(textscan(fid, '%f')).';
        fclose(fid);  
        hcval_path_sf = reshape(hcval_path_sf,SS,JJ+1);

   fid = fopen('hc_return\simdata\hcval_path_sm.txt', 'r');
        hcval_path_sm  = cell2mat(textscan(fid, '%f')).';
        fclose(fid);  
        hcval_path_sm = reshape(hcval_path_sm,SS,JJ+1);

    % import simulated risky share (single women and single men)
    fid = fopen('hc_return\simdata\share_path_sf.txt', 'r');
        share_path_sf  = cell2mat(textscan(fid, '%f')).';
        fclose(fid);  
        share_path_sf = reshape(share_path_sf,SS,JJ+1);

   fid = fopen('hc_return\simdata\share_path_sm.txt', 'r');
        share_path_sm  = cell2mat(textscan(fid, '%f')).';
        fclose(fid);  
        share_path_sm = reshape(share_path_sm,SS,JJ+1);
  
   % transform negative values to NaN (refers to observations that are not single)     
   hcval_path_sm(hcval_path_sm < 0) = nan;
   hcval_path_sf(hcval_path_sf < 0) = nan;

   % Compute income for single men and single women

   % Import income components: single women 
    fid = fopen('hc_return\yfs.txt', 'r');
        yfs  = cell2mat(textscan(fid, '%f')).';
        fclose(fid);  

    fid = fopen('hc_return\efffs.txt', 'r');
        efffs  = cell2mat(textscan(fid, '%f')).';
        fclose(fid); 
        
    fid = fopen('hc_return\thetafs.txt', 'r');
        thetafs  = cell2mat(textscan(fid, '%f')).';
        fclose(fid); 

    fid = fopen('hc_return\afs_boom.txt', 'r');
        afs_boom  = cell2mat(textscan(fid, '%f')).';
        fclose(fid); 

    fid = fopen('hc_return\afs_rec.txt', 'r');
        afs_rec  = cell2mat(textscan(fid, '%f')).';
        fclose(fid); 

    % Import income components: single men  
    fid = fopen('hc_return\yms.txt', 'r');
        yms  = cell2mat(textscan(fid, '%f')).';
        fclose(fid);  

    fid = fopen('hc_return\effms.txt', 'r');
        effms  = cell2mat(textscan(fid, '%f')).';
        fclose(fid); 
        
    fid = fopen('hc_return\thetams.txt', 'r');
        thetams  = cell2mat(textscan(fid, '%f')).';
        fclose(fid); 

    fid = fopen('hc_return\ams_boom.txt', 'r');
        ams_boom  = cell2mat(textscan(fid, '%f')).';
        fclose(fid); 

    fid = fopen('hc_return\ams_rec.txt', 'r');
        ams_rec  = cell2mat(textscan(fid, '%f')).';
        fclose(fid); 

    % compute chain for stochastic part of income 
    % single women
    chain_sincf = zeros(SS,JJ);
    for t = 1: JJ+1
    for s = 1:SS
    if (chain_agg(s,t) == 1) 
        chain_sincf(s,t) = afs_boom(chain_prodf(s,t));
    elseif (chain_agg(s,t) == 2) 
    chain_sincf(s,t) = afs_rec(chain_prodf(s,t));
    end 
    end  
    end  
    % set couple observations to NaN
    chain_sincf(chain_martf>0) = NaN;

    % single men
    chain_sincm = zeros(SS,JJ);
    for t = 1: JJ+1
    for s = 1:SS
    if (chain_agg(s,t) == 1) 
        chain_sincm(s,t) = ams_boom(chain_prodm(s,t));
    elseif (chain_agg(s,t) == 2) 
    chain_sincm(s,t) = ams_rec(chain_prodm(s,t));
    end 
    end  
    end  
    % set couple observations to NaN
    chain_sincm(chain_martm>0) = NaN;

    % compute pathsf of human capital returns (single men and single women)
    hcret_path_sm = NaN(SS,JJ);
    hcret_path_sf = NaN(SS,JJ);

    % share of high and low education 
    AB_f = round(0.44*SS);
    AB_m = round(0.41*SS);

  for t = 2: JJ
   
    for s = 1:AB_f
     hcret_path_sf(s,t) = (hcval_path_sf(s,t) + yfs*efffs(t)*thetafs(1)*chain_sincf(s,t))/hcval_path_sf(s,t-1);
    end
    for s = AB_f+1:SS
     hcret_path_sf(s,t) = (hcval_path_sf(s,t) + yfs*efffs(t)*thetafs(2)*chain_sincf(s,t))/hcval_path_sf(s,t-1);
    end
    for s = 1:AB_m
     hcret_path_sm(s,t) = (hcval_path_sm(s,t) + yms*effms(t)*thetams(1)*chain_sincm(s,t))/hcval_path_sm(s,t-1);
    end
    for s = AB_m+1:SS
     hcret_path_sm(s,t) = (hcval_path_sm(s,t) + yms*effms(t)*thetams(2)*chain_sincm(s,t))/hcval_path_sm(s,t-1);
    end

  end

    % compute path of realized stock returns (single men and single women)
    fid = fopen('hc_return\ard_boom.txt', 'r');
        ard_boom  = cell2mat(textscan(fid, '%f')).';
        fclose(fid); 

    fid = fopen('hc_return\ard_rec.txt', 'r');
        ard_rec  = cell2mat(textscan(fid, '%f')).';
        fclose(fid); 

    ret_realf = zeros(SS,JJ);
    ret_realm = zeros(SS,JJ);

    for t = 1: JJ
        for s = 1:SS
   if(chain_agg(s,t) == 1)
   ret_realf(s,t) = ard_boom(chain_retf(s,t));
   ret_realm(s,t) = ard_boom(chain_retm(s,t));
   elseif(chain_agg(s,t) == 2)
   ret_realf(s,t) = ard_rec(chain_retf(s,t));
   ret_realm(s,t) = ard_rec(chain_retm(s,t));
   end 
        end
    end

  % set couple observations to NaN
   ret_realf(isnan(hcret_path_sf))= nan;
   ret_realm(isnan(hcret_path_sm))= nan;

  % compute correlation of human capital returns and stock market returns
   corr_hcretf = NaN(1,JJ);
   corr_hcretm = NaN(1,JJ);

    for t = 1: JJ
       corr_hcretf(t) = corr(hcret_path_sf(:,t),ret_realf(:,t),'rows','complete');
       corr_hcretm(t) = corr(hcret_path_sm(:,t),ret_realm(:,t),'rows','complete');
    end   

    % smooth figures (moving average of 3)
    corr_hcretf_smooth = movmean(corr_hcretf,3);
    corr_hcretm_smooth = movmean(corr_hcretm,3);

    % plot results: figure 13 in paper
        figure(13)
        hold on
        plot(age,   corr_hcretf_smooth, 'Linewidth',  3.5,'LineStyle', '--','Color','black')
        plot(age,   corr_hcretm_smooth, 'Linewidth',  3.5,'Color','black')
        hold off
        xlim([32 65])
        ylim([-0.1 0.3])
        xlabel('Age','Interpreter','LaTex','Fontsize',24)
        set(gca,'XGrid','off','YGrid','on','Fontsize',24)
        leg = legend('single women','single men','Interpreter','LaTex');
        set(leg,'Interpreter','LaTex','Fontsize',24,'Location','NorthEast')
        saveas(gcf,[figurefolder '\hcret_correlation.png'])
        hold off

%% prepare counterfactuals 1e and 1f below
  
% Deterministic income components: single women 
    fid = fopen('base\yfs.txt', 'r');
        yfs  = cell2mat(textscan(fid, '%f')).';
        fclose(fid);  

    fid = fopen('base\efffs.txt', 'r');
        efffs  = cell2mat(textscan(fid, '%f')).';
        fclose(fid); 
        
    fid = fopen('base\thetafs.txt', 'r');
        thetafs  = cell2mat(textscan(fid, '%f')).';
        fclose(fid); 

    % Deterministic income components: single men  
    fid = fopen('base\yms.txt', 'r');
        yms  = cell2mat(textscan(fid, '%f')).';
        fclose(fid);  

    fid = fopen('base\effms.txt', 'r');
        effms  = cell2mat(textscan(fid, '%f')).';
        fclose(fid); 
        
    fid = fopen('base\thetams.txt', 'r');
        thetams  = cell2mat(textscan(fid, '%f')).';
        fclose(fid); 

 % share of high and low education 
 AB_f = round(0.44*SS);
 AB_m = round(0.41*SS);

% compute deterministic income part for single men and single women
det_inc_sf = zeros(SS,JJ);
det_inc_sm = zeros(SS,JJ);

for t = 1: JJ 
    for s = 1:AB_f
    det_inc_sf(s,t) = yfs*efffs(t)*thetafs(1);
    end
    for s = AB_f+1:SS
    det_inc_sf(s,t) = yfs*efffs(t)*thetafs(2);
    end
    for s = 1:AB_m
    det_inc_sm(s,t) = yms*effms(t)*thetams(1);
    end
    for s = AB_m+1:SS
    det_inc_sm(s,t) = yms*effms(t)*thetams(2);
    end
end

% COUNTERFACTUAL 1e
% if I change age term, what is the correction needed to keep same level as
% in baseline?

% idea behind this counterfactual: fix the female deterministic income level, but change their slope to male
% That is: on average, women have the same deterministic income as before.

% counterfactual female income: female level and education term, but male
% age term
det_inc_sf_cnt = zeros(SS,JJ);

for t = 1: JJ 
    for s = 1:AB_f
    det_inc_sf_cnt(s,t) = yfs*effms(t)*thetafs(1);
    end
    for s = AB_f+1:SS
    det_inc_sf_cnt(s,t) = yfs*effms(t)*thetafs(2);
    end
end

% what is the average difference between actual and counterfactual det.
% income?
adj_term = mean(mean(det_inc_sf))/mean(mean(det_inc_sf_cnt));

% compute counterfactual female income path that assigns women male slopes
% but then corrects level "os: only slope"
det_inc_sf_os = zeros(SS,JJ);

for t = 1: JJ 
    for s = 1:AB_f
    det_inc_sf_os(s,t) = yfs*effms(t)*thetafs(1)*adj_term ;
    end
    for s = AB_f+1:SS
    det_inc_sf_os(s,t) = yfs*effms(t)*thetafs(2)*adj_term ;
    end
end

% check
mean(mean(det_inc_sf))
mean(mean(det_inc_sf_cnt))
mean(mean(det_inc_sf_os))

% COUNTERFACTUAL 1f
% if I change constant term, what is the correction needed to keep slope as in
% baseline
% idea behind this counterfactual: adjust female det. income level to be same as for
% men, but keep the female slope
% That is: on average, men and women have the same deterministic income now.  

% what is the avaerge difference between male and female det. income?
diff_inc = mean(mean(det_inc_sm))-mean(mean(det_inc_sf));

% compute counterfactual income profile "ol = only level"
% by shifting female profile up to be on average same as men
det_inc_sf_ol = zeros(SS,JJ);
for t = 1: JJ 
    for s = 1:AB_f
    det_inc_sf_ol(s,t) = (yfs*efffs(t)*thetafs(1))+diff_inc;
    end
    for s = AB_f+1:SS
    det_inc_sf_ol(s,t) = (yfs*efffs(t)*thetafs(2))+diff_inc;
    end
end

% check
mean(mean(det_inc_sm))
mean(mean(det_inc_sf_ol))

% Note: cannot compute like that in fortran code 
% because I work with individual elements, and they are multiplicative
% therefore: check that I get the same results when done as in fortran code (multiplicative instead of
% additive)

det_inc_sf_ol_fort = zeros(SS,JJ);
add_low=zeros(JJ);
add_high=zeros(JJ);
efffs_ol= zeros(JJ);

% relative increase for high and low educated 
for j=1:JJ
add_low(j) = ((yfs*efffs(j)*thetafs(1))+diff_inc)/(yfs*efffs(j)*thetafs(1));
add_high(j)= ((yfs*efffs(j)*thetafs(2))+diff_inc)/(yfs*efffs(j)*thetafs(2));
end

% adjust age term by relative income for low
for j=1:JJ
efffs_ol(j) = efffs(j)*add_low(j);
end 

% then, adjust high education term to correct for differences across
% education levels 
% Note: this correction should be age-dependent, but when doing it for age 10:
% approximation is very close (see graph below)
thetafs_ol(2) = thetafs(2) *(add_high(10)/add_low(10));

% counterfactual determinstic income path
for t = 1: JJ 
    for s = 1:AB_f
    det_inc_sf_ol_fort(s,t) = (yfs*efffs_ol(t)*thetafs(1));
    end
    for s = AB_f+1:SS
    det_inc_sf_ol_fort(s,t) = (yfs*efffs_ol(t)*thetafs_ol(2));
    end
end

% Check for approximation: very close when doing it for age 10
hold on
plot(age, mean(det_inc_sf_ol_fort), 'Linewidth',  3.5,'LineStyle', '--','Color','black')
plot(age, mean(det_inc_sf_ol), 'Linewidth',  3.5,'Color','black')
hold off

        
%% plot data model comparison of life-cycle profiles (Figures 3,4, and 14)

        % Stock Market Participation Rate
        % Figure 3b
        figure(32)  
        hold on
        plot(age,  SMPcoh_sf,'Linewidth', 3.5,'LineStyle', '--','Color','black')
        plot(age,  SMPcoh_sf_dat, 'Linewidth', 1.5, 'Color','black')
        plot(age,  SMPcoh_sf_ch_dat,'Linewidth', .5,'color','.5 .5 .5')
        plot(age,  SMPcoh_sf_cl_dat,'Linewidth', .5,'Color','.5 .5 .5')
        hold off
        xlim([30 64])
        ylim([0 1])
        xlabel('Age','Interpreter','LaTex','Fontsize',32)
        ylabel('Participation Rate','Interpreter','LaTex','Fontsize',32)
        set(gca,'XGrid','off','YGrid','on','Fontsize',24)
        saveas(gcf,[figurefolder '\SMP_lifecycle_model_sf.png'])
        hold off
        
        % Figure 3e
        figure(35)
        hold on
        plot(age,  SMPcoh_sm_dat, 'Linewidth', 1.5, 'Color','black')
        plot(age,  SMPcoh_sm_ch_dat,'Linewidth', .5,'color','.5 .5 .5')
        plot(age,  SMPcoh_sm_cl_dat,'Linewidth', .5,'Color','.5 .5 .5')
        plot(age,  SMPcoh_sm,'Linewidth', 3.5,'LineStyle', '--','Color','black')
        hold off
        xlim([30 64])
        ylim([0 1])
        xlabel('Age','Interpreter','LaTex','Fontsize',32)
        ylabel('Participation Rate','Interpreter','LaTex','Fontsize',32)
        set(gca,'XGrid','off','YGrid','on','Fontsize',24)
        saveas(gcf,[figurefolder '\SMP_lifecycle_model_sm.png'])
        hold off
        
        % Figure 14b
        figure(142)
        hold on
        plot(age,  SMPcoh_c ,'Linewidth', 3.5, 'LineStyle', '--','Color','black')
        plot(age,  SMP_c_dat, 'Linewidth', 1.5,'Color','black')
        plot(age,  SMPcoh_c_ch_dat,'Linewidth', .5,'color','.5 .5 .5')
        plot(age,  SMPcoh_c_cl_dat,'Linewidth', .5,'Color','.5 .5 .5')
        hold off
        xlim([30 64])
        ylim([0 1])
        xlabel('Age','Interpreter','LaTex','Fontsize',32)
        ylabel('Participation Rate','Interpreter','LaTex','Fontsize',32)
        leg = legend('model','data','Interpreter','LaTex');
        set(leg,'Interpreter','LaTex','Fontsize',30,'Location','NorthEast','NumColumns',2);
        set(gca,'XGrid','off','YGrid','on','Fontsize',24)
        saveas(gcf,[figurefolder '\SMP_lifecycle_model_c.png'])
        hold off
        
        % Net Worth 
        % Figure 4a
        figure(41)
        hold on
        plot(age,  cashcoh_sf*5, 'Linewidth', 3.5, 'LineStyle', '--','Color','black')
        plot(age,  nw_sf_dat_smooth,'Linewidth', 1.5,'Color','black')
        plot(age,  nw_sf_ch_sm,'Linewidth', .5,'color','.5 .5 .5')
        plot(age,  nw_sf_cl_sm,'Linewidth', .5,'Color','.5 .5 .5')
        hold off
        xlim([30 64])
        ylim([0 600])
        yticks(0:200:600)
        xlabel('Age','Interpreter','LaTex','Fontsize',32)
        ylabel('(000s) Assets in 2007 \$','Interpreter','LaTex','Fontsize',32)
        set(gca,'XGrid','off','YGrid','on','Fontsize',24)
        saveas(gcf,[figurefolder '\nw_lifecycle_model_sf.png'])
        hold off
        
        % Figure 4b
        figure(42)
        hold on
        plot(age,  cashcoh_sm*5, 'Linewidth', 3.5, 'LineStyle', '--','Color','black')
        plot(age,  nw_sm_dat_smooth,'Linewidth', 1.5,'Color','black')
        plot(age,  nw_sm_ch_sm,'Linewidth', .5,'color','.5 .5 .5')
        plot(age,  nw_sm_cl_sm,'Linewidth', .5,'Color','.5 .5 .5')
        hold off
        xlim([30 64])
        ylim([0 600])
        yticks(0:200:600)
        xlabel('Age','Interpreter','LaTex','Fontsize',32)
        ylabel('(000s) Assets in 2007 \$','Interpreter','LaTex','Fontsize',32)
        leg = legend('model','data','Interpreter','LaTex');
        set(leg,'Interpreter','LaTex','Fontsize',30,'Location','NorthEast','NumColumns',2)
        set(gca,'XGrid','off','YGrid','on','Fontsize',24)
        saveas(gcf,[figurefolder '\nw_lifecycle_model_sm.png'])
        hold off
        
        % Figure 14d
        figure(144)
        hold on
        plot(age,  cashcoh_c*5,'Linewidth', 3.5,'LineStyle', '--','Color','black')
        plot(age,  nw_c_dat, 'Linewidth', 1.5,'Color', 'black')
        plot(age,  nw_c_ch_dat,'Linewidth', .5,'color','.5 .5 .5')
        plot(age,  nw_c_cl_dat,'Linewidth', .5,'Color','.5 .5 .5')
        hold off
        xlim([30 64])
        ylim([0 1000])
        xlabel('Age','Interpreter','LaTex','Fontsize',32)
        ylabel('(000s) Assets in 2007 \$','Interpreter','LaTex','Fontsize',32)
        set(gca,'XGrid','off','YGrid','on','Fontsize',24)
        saveas(gcf,[figurefolder '\nw_lifecycle_model_c.png'])
        hold off
    
        % Conditional Risky Share
        % Figure 3f
        figure(36)
        hold on
        plot(age,  sharecoh_sm,'Linewidth', 3.5,'Color','black','LineStyle', '--')
        plot(age,  sharecoh_sm_dat, 'Linewidth', 1.5, 'Color','black')
        plot(age,  sharecoh_sm_ch_dat,'Linewidth', .5,'color','.5 .5 .5')
        plot(age,  sharecoh_sm_cl_dat,'Linewidth', .5,'Color','.5 .5 .5')
        hold off
        xlim([30 64])
        ylim([0 1])
        xlabel('Age','Interpreter','LaTex','Fontsize',24)
        ylabel('Equity Share','Interpreter','LaTex','Fontsize',32)
        set(gca,'XGrid','off','YGrid','on','Fontsize',24)
        saveas(gcf,[figurefolder '\share_lifecycle_model_sm.png'])
        hold off
        
        % Figure 3c
        figure(33)
        hold on
        plot(age,  sharecoh_sf,'Linewidth', 3.5,'Color','black','LineStyle', '--')
        plot(age,  sharecoh_sf_dat, 'Linewidth', 1.5, 'Color','black')
        plot(age,  sharecoh_sf_ch_dat,'Linewidth', .5,'color','.5 .5 .5')
        plot(age,  sharecoh_sf_cl_dat,'Linewidth', .5,'Color','.5 .5 .5')
        hold off
        xlim([30 64])
        ylim([0 1])
        xlabel('Age','Interpreter','LaTex','Fontsize',24)
        ylabel('Equity Share','Interpreter','LaTex','Fontsize',32)
        set(gca,'XGrid','off','YGrid','on','Fontsize',24)
        leg = legend('model','data','Interpreter','LaTex');
        set(leg,'Interpreter','LaTex','Fontsize',32,'Location','NorthEast','NumColumns',2);
        saveas(gcf,[figurefolder '\share_lifecycle_model_sf.png'])
        hold off

        % Figure 14c
        figure(143)
        hold on
        plot(age,  sharecoh_c,'Linewidth', 3.5,'Color','black','LineStyle', '--')
        plot(age,  sharecoh_c_dat, 'Linewidth', 1.5, 'Color','black')
        plot(age,  sharecoh_c_ch_dat,'Linewidth', .5,'color','.5 .5 .5')
        plot(age,  sharecoh_c_cl_dat,'Linewidth', .5,'Color','.5 .5 .5')
        hold off
        xlim([30 64])
        ylim([0 1])
        xlabel('Age','Interpreter','LaTex','Fontsize',24)
        ylabel('Equity Share','Interpreter','LaTex','Fontsize',32)
        set(gca,'XGrid','off','YGrid','on','Fontsize',24)
        saveas(gcf,[figurefolder '\share_lifecycle_model_c.png'])
        hold off
        
        % Unconditional Share
        % Figure 3d
        figure(34)
        hold on
        plot(age,  ucsharecoh_sm,'Linewidth', 3.5,'Color','black','LineStyle', '--')
        plot(age,  ucsharecoh_sm_dat, 'Linewidth', 1.5, 'Color','black')
        plot(age,  ucsharecoh_sm_ch_dat,'Linewidth', .5,'color','.5 .5 .5')
        plot(age,  ucsharecoh_sm_cl_dat,'Linewidth', .5,'Color','.5 .5 .5')
        hold off
        xlim([30 64])
        ylim([0 0.5])
        yticks(0:0.1:0.5);
        xlabel('Age','Interpreter','LaTex','Fontsize',24)
        ylabel('Equity Share','Interpreter','LaTex','Fontsize',32)
        set(gca,'XGrid','off','YGrid','on','Fontsize',24)
        saveas(gcf,[figurefolder '\ucshare_lifecycle_model_sm.png'])
        hold off

        % Figure 3a
        figure(31)
        hold on
        plot(age,  ucsharecoh_sf,'Linewidth', 3.5,'Color','black','LineStyle', '--')
        plot(age,  ucsharecoh_sf_dat, 'Linewidth', 1.5, 'Color','black')
        plot(age,  ucsharecoh_sf_ch_dat,'Linewidth', .5,'color','.5 .5 .5')
        plot(age,  ucsharecoh_sf_cl_dat,'Linewidth', .5,'Color','.5 .5 .5')
        hold off
        xlim([30 64])
        ylim([0 0.5])
        yticks(0:0.1:0.5)
        xlabel('Age','Interpreter','LaTex','Fontsize',24)
        ylabel('Equity Share','Interpreter','LaTex','Fontsize',32)
        set(gca,'XGrid','off','YGrid','on','Fontsize',24)
        saveas(gcf,[figurefolder '\ucshare_lifecycle_model_sf.png'])
        hold off

        % Figure 14a
        figure(141)
        hold on
        plot(age,  ucsharecoh_c,'Linewidth', 3.5,'Color','black','LineStyle', '--')
        plot(age,  ucsharecoh_c_dat, 'Linewidth', 1.5, 'Color','black')
        plot(age,  ucsharecoh_c_ch_dat,'Linewidth', .5,'color','.5 .5 .5')
        plot(age,  ucsharecoh_c_cl_dat,'Linewidth', .5,'Color','.5 .5 .5')
        hold off
        xlim([30 64])
        ylim([0 0.5])
        yticks(0:0.1:0.5);
        xlabel('Age','Interpreter','LaTex','Fontsize',24)
        ylabel('Equity Share','Interpreter','LaTex','Fontsize',32)
        set(gca,'XGrid','off','YGrid','on','Fontsize',24)
        saveas(gcf,[figurefolder '\ucshare_lifecycle_model_c.png'])
        hold off


        %% %%%%%%%%%%%%% decomposition (Section 6.1) %%%%%%%%%%%%%%%%%%

%% decomposition 1:  import results deterministic income
   
    % EQUITY SHARE
    fid = fopen('decomp1\cohort\ucsharecoh_sf.txt');
    ucshare_sf_dc1  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('decomp1\cohort\ucsharecoh_sm.txt');
    ucshare_sm_dc1  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('decomp1\cohort\ucsharecoh_c.txt');
    ucshare_c_dc1   = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    % ASSET HOLDINGS
    fid = fopen('decomp1\cohort\cashcoh_sf.txt');
    cash_sf_dc1  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('decomp1\cohort\cashcoh_sm.txt');
    cash_sm_dc1  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('decomp1\cohort\cashcoh_c.txt');
    cash_c_dc1   = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

%% decomposition 1b: import results det. income: only fixed term

    % EQUITY SHARE
    fid = fopen('decomp1b\cohort\ucsharecoh_sf.txt');
    ucshare_sf_dc1b  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('decomp1b\cohort\ucsharecoh_sm.txt');
    ucshare_sm_dc1b  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('decomp1b\cohort\ucsharecoh_c.txt');
    ucshare_c_dc1b   = cell2mat(textscan(fid, '%f')).';
    fclose(fid);
    
    % ASSET HOLDINGS
    fid = fopen('decomp1b\cohort\cashcoh_sf.txt');
    cash_sf_dc1b  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('decomp1b\cohort\cashcoh_sm.txt');
    cash_sm_dc1b  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('decomp1b\cohort\cashcoh_c.txt');
    cash_c_dc1b   = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

%% decomposition 1c: import results det. income: only age term

    % EQUITY SHARE
    fid = fopen('decomp1c\cohort\ucsharecoh_sf.txt');
    ucshare_sf_dc1c  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('decomp1c\cohort\ucsharecoh_sm.txt');
    ucshare_sm_dc1c  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('decomp1c\cohort\ucsharecoh_c.txt');
    ucshare_c_dc1c   = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    % ASSET HOLDINGS
    fid = fopen('decomp1c\cohort\cashcoh_sf.txt');
    cash_sf_dc1c  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('decomp1c\cohort\cashcoh_sm.txt');
    cash_sm_dc1c  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('decomp1c\cohort\cashcoh_c.txt');
    cash_c_dc1c   = cell2mat(textscan(fid, '%f')).';
    fclose(fid);      

%% decomposition 1d: import results det. income: only educ term

    % EQUITY SHARE
    fid = fopen('decomp1d\cohort\ucsharecoh_sf.txt');
    ucshare_sf_dc1d  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('decomp1d\cohort\ucsharecoh_sm.txt');
    ucshare_sm_dc1d  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('decomp1d\cohort\ucsharecoh_c.txt');
    ucshare_c_dc1d   = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    % ASSET HOLDINGS
    fid = fopen('decomp1d\cohort\cashcoh_sf.txt');
    cash_sf_dc1d  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('decomp1d\cohort\cashcoh_sm.txt');
    cash_sm_dc1d  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('decomp1d\cohort\cashcoh_c.txt');
    cash_c_dc1d   = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 

%% decomposition 1e: import results det. income: same level, different slope

    % EQUITY SHARE
    fid = fopen('decomp1e\cohort\ucsharecoh_sf.txt');
    ucshare_sf_dc1e  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('decomp1e\cohort\ucsharecoh_sm.txt');
    ucshare_sm_dc1e  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('decomp1e\cohort\ucsharecoh_c.txt');
    ucshare_c_dc1e   = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    % ASSET HOLDINGS
    fid = fopen('decomp1e\cohort\cashcoh_sf.txt');
    cash_sf_dc1e  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('decomp1e\cohort\cashcoh_sm.txt');
    cash_sm_dc1e  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('decomp1e\cohort\cashcoh_c.txt');
    cash_c_dc1e   = cell2mat(textscan(fid, '%f')).';
    fclose(fid);  

%% decomposition 1f: import results det. income: same slope, different level

    % EQUITY SHARE
    fid = fopen('decomp1f\cohort\ucsharecoh_sf.txt');
    ucshare_sf_dc1f  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('decomp1f\cohort\ucsharecoh_sm.txt');
    ucshare_sm_dc1f  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('decomp1f\cohort\ucsharecoh_c.txt');
    ucshare_c_dc1f   = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    % ASSET HOLDINGS
    fid = fopen('decomp1f\cohort\cashcoh_sf.txt');
    cash_sf_dc1f  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('decomp1f\cohort\cashcoh_sm.txt');
    cash_sm_dc1f  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('decomp1f\cohort\cashcoh_c.txt');
    cash_c_dc1f   = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 

%% decomposition 2:  import results stochastic income 

    % EQUITY SHARE
    fid = fopen('decomp2\cohort\ucsharecoh_sf.txt');
    ucshare_sf_dc2  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('decomp2\cohort\ucsharecoh_sm.txt');
    ucshare_sm_dc2  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('decomp2\cohort\ucsharecoh_c.txt');
    ucshare_c_dc2   = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    % ASSET HOLDINGS
    fid = fopen('decomp2\cohort\cashcoh_sf.txt');
    cash_sf_dc2  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('decomp2\cohort\cashcoh_sm.txt');
    cash_sm_dc2  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('decomp2\cohort\cashcoh_c.txt');
    cash_c_dc2   = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

%% decomposition 3:  import results household sizes 
    
    % EQUITY SHARE

    fid = fopen('decomp3\cohort\ucsharecoh_sf.txt');
    ucshare_sf_dc3  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('decomp3\cohort\ucsharecoh_sm.txt');
    ucshare_sm_dc3  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('decomp3\cohort\ucsharecoh_c.txt');
    ucshare_c_dc3  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    % ASSET HOLDINGS
    fid = fopen('decomp3\cohort\cashcoh_sf.txt');
    cash_sf_dc3  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('decomp3\cohort\cashcoh_sm.txt');
    cash_sm_dc3  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('decomp3\cohort\cashcoh_c.txt');
    cash_c_dc3   = cell2mat(textscan(fid, '%f')).';
    fclose(fid);
 
%% decomposition 4:  import results marriage probabilities
    
    % EQUITY SHARE
    fid = fopen('decomp4\cohort\ucsharecoh_sf.txt');
    ucshare_sf_dc4  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('decomp4\cohort\ucsharecoh_sm.txt');
    ucshare_sm_dc4  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('decomp4\cohort\ucsharecoh_c.txt');
    ucshare_c_dc4  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    % ASSET HOLDINGS
    fid = fopen('decomp4\cohort\cashcoh_sf.txt');
    cash_sf_dc4  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('decomp4\cohort\cashcoh_sm.txt');
    cash_sm_dc4  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('decomp4\cohort\cashcoh_c.txt');
    cash_c_dc4   = cell2mat(textscan(fid, '%f')).';
    fclose(fid);
        
%% decomposition 5:  import results marriage mapping 
    
    % EQUITY SHARE
    fid = fopen('decomp5\cohort\ucsharecoh_sf.txt');
    ucshare_sf_dc5  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('decomp5\cohort\ucsharecoh_sm.txt');
    ucshare_sm_dc5  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('decomp5\cohort\ucsharecoh_c.txt');
    ucshare_c_dc5  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    % ASSET HOLDINGS
    fid = fopen('decomp5\cohort\cashcoh_sf.txt');
    cash_sf_dc5  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('decomp5\cohort\cashcoh_sm.txt');
    cash_sm_dc5  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('decomp5\cohort\cashcoh_c.txt');
    cash_c_dc5   = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

%% decomposition 6:  import results initial wealth holdings
    
    % EQUITY SHARE
    fid = fopen('decomp6\cohort\ucsharecoh_sf.txt');
    ucshare_sf_dc6  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('decomp6\cohort\ucsharecoh_sm.txt');
    ucshare_sm_dc6  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('decomp6\cohort\ucsharecoh_c.txt');
    ucshare_c_dc6  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    % ASSET HOLDINGS
    fid = fopen('decomp6\cohort\cashcoh_sf.txt');
    cash_sf_dc6  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('decomp6\cohort\cashcoh_sm.txt');
    cash_sm_dc6  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('decomp6\cohort\cashcoh_c.txt');
    cash_c_dc6   = cell2mat(textscan(fid, '%f')).';
    fclose(fid);
      
%% decomposition 7:  import results education distribution
    
    % EQUITY SHARE
    fid = fopen('decomp7\cohort\ucsharecoh_sf.txt');
    ucshare_sf_dc7  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('decomp7\cohort\ucsharecoh_sm.txt');
    ucshare_sm_dc7  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('decomp7\cohort\ucsharecoh_c.txt');
    ucshare_c_dc7  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    % ASSET HOLDINGS
    fid = fopen('decomp7\cohort\cashcoh_sf.txt');
    cash_sf_dc7  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('decomp7\cohort\cashcoh_sm.txt');
    cash_sm_dc7  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('decomp7\cohort\cashcoh_c.txt');
    cash_c_dc7   = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

%% decomposition 8:  import results medical expenses

    % EQUITY SHARE
    fid = fopen('decomp8\cohort\ucsharecoh_sf.txt');
    ucshare_sf_dc8  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('decomp8\cohort\ucsharecoh_sm.txt');
    ucshare_sm_dc8  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('decomp8\cohort\ucsharecoh_c.txt');
    ucshare_c_dc8  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    % ASSET HOLDINGS
    fid = fopen('decomp8\cohort\cashcoh_sf.txt');
    cash_sf_dc8  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('decomp8\cohort\cashcoh_sm.txt');
    cash_sm_dc8  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('decomp8\cohort\cashcoh_c.txt');
    cash_c_dc8   = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

%% decomposition 9:  import results survival risk

    % EQUITY SHARE
    fid = fopen('decomp9\cohort\ucsharecoh_sf.txt');
    ucshare_sf_dc9  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('decomp9\cohort\ucsharecoh_sm.txt');
    ucshare_sm_dc9  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('decomp9\cohort\ucsharecoh_c.txt');
    ucshare_c_dc9  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    % ASSET HOLDINGS
    fid = fopen('decomp9\cohort\cashcoh_sf.txt');
    cash_sf_dc9  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('decomp9\cohort\cashcoh_sm.txt');
    cash_sm_dc9  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('decomp9\cohort\cashcoh_c.txt');
    cash_c_dc9   = cell2mat(textscan(fid, '%f')).';
    fclose(fid);


%% combine all decomposition results: Table 5 in paper

% compute gaps in data (equity shares and asset holdings)
 ucshare_gap_dat = mean(ucsharecoh_sm_dat) - mean(ucsharecoh_sf_dat); 
 cash_gap_dat    = mean(nw_sm_dat_smooth) - mean(nw_sf_dat_smooth);

 % compute gap in asset holdings in model: baseline and all counterfactuals
    cash_gap = NaN(1,15);
    cash_gap(1)  = (mean(cashcoh_sm)  -  mean(cashcoh_sf))*5; 
    cash_gap(2)  = (mean(cash_sm_dc1) -  mean(cash_sf_dc1))*5;
    cash_gap(3)  = (mean(cash_sm_dc1b)-  mean(cash_sf_dc1b))*5;
    cash_gap(4)  = (mean(cash_sm_dc1c)-  mean(cash_sf_dc1c))*5;
    cash_gap(5)  = (mean(cash_sm_dc1d)-  mean(cash_sf_dc1d))*5;
    cash_gap(6)  = (mean(cash_sm_dc1e)-  mean(cash_sf_dc1e))*5;
    cash_gap(7)  = (mean(cash_sm_dc1f)-  mean(cash_sf_dc1f))*5;
    cash_gap(8)  = (mean(cash_sm_dc2) -  mean(cash_sf_dc2))*5;  
    cash_gap(9)  = (mean(cash_sm_dc3) -  mean(cash_sf_dc3))*5;
    cash_gap(10) = (mean(cash_sm_dc4) -  mean(cash_sf_dc4))*5;
    cash_gap(11) = (mean(cash_sm_dc5) -  mean(cash_sf_dc5))*5;
    cash_gap(12) = (mean(cash_sm_dc6) -  mean(cash_sf_dc6))*5;
    cash_gap(13) = (mean(cash_sm_dc7) -  mean(cash_sf_dc7))*5;
    cash_gap(14) = (mean(cash_sm_dc8) -  mean(cash_sf_dc8))*5;
    cash_gap(15) = (mean(cash_sm_dc9) -  mean(cash_sf_dc9))*5;
    
    X= ['Decomp asset:   ','baseline: ' , num2str(cash_gap(1)),' income level: ' , ...
           num2str(cash_gap(2)), '  fixed term: ' , num2str(cash_gap(3)), ' age term: ' , num2str(cash_gap(4)), ...
         ' educ term: ' , num2str(cash_gap(5)),' same level (slope effect): ' , num2str(cash_gap(6)), ...
         ' same slope (level effect): ' , num2str(cash_gap(7)),' income risk: ' , num2str(cash_gap(8)),...
         ' HH size: ' , num2str(cash_gap(9)),' marriage prob: ' , num2str(cash_gap(10)),... 
         ' marriage market: ' , num2str(cash_gap(11)), ' initial wealth: ' , num2str(cash_gap(12)),...
         ' education distribution: ' , num2str(cash_gap(13)),' medical expenses: ' , num2str(cash_gap(14)),...
         ' survival risk: ' , num2str(cash_gap(15))];
    disp(X);

    % what part of gap can be explained by individual channels?
    cash_gap_expl = NaN(1,15);
    cash_gap_expl(1)  = 1-(cash_gap(1)/cash_gap(1));
    cash_gap_expl(2)  = 1-(cash_gap(2)/cash_gap(1));
    cash_gap_expl(3)  = 1-(cash_gap(3)/cash_gap(1));
    cash_gap_expl(4)  = 1-(cash_gap(4)/cash_gap(1));
    cash_gap_expl(5)  = 1-(cash_gap(5)/cash_gap(1));
    cash_gap_expl(6)  = 1-(cash_gap(6)/cash_gap(1));
    cash_gap_expl(7)  = 1-(cash_gap(7)/cash_gap(1));
    cash_gap_expl(8)  = 1-(cash_gap(8)/cash_gap(1));
    cash_gap_expl(9)  = 1-(cash_gap(9)/cash_gap(1));
    cash_gap_expl(10) = 1-(cash_gap(10)/cash_gap(1));
    cash_gap_expl(11) = 1-(cash_gap(11)/cash_gap(1));
    cash_gap_expl(12) = 1-(cash_gap(12)/cash_gap(1));
    cash_gap_expl(13) = 1-(cash_gap(13)/cash_gap(1));
    cash_gap_expl(14) = 1-(cash_gap(14)/cash_gap(1));
    cash_gap_expl(15) = 1-(cash_gap(15)/cash_gap(1));
     
  X= ['Decomp asset explained:   ','baseline: ' , num2str(cash_gap_expl(1)),' income level: ' , ...
           num2str(cash_gap_expl(2)), '  fixed term: ' , num2str(cash_gap_expl(3)), ' age term: ' , num2str(cash_gap_expl(4)), ...
         ' educ term: ' , num2str(cash_gap_expl(5)),' only level: ' , num2str(cash_gap_expl(6)), ...
         ' only slope: ' , num2str(cash_gap_expl(7)),' income risk: ' , num2str(cash_gap_expl(8)),...
         ' HH size: ' , num2str(cash_gap_expl(9)),' marriage prob: ' , num2str(cash_gap_expl(10)),... 
         ' marriage market: ' , num2str(cash_gap_expl(11)), ' initial wealth: ' , num2str(cash_gap_expl(12)),...
         ' education distribution: ' , num2str(cash_gap_expl(13)),' medical expenses: ' , num2str(cash_gap_expl(14)),...
         ' survival risk: ' , num2str(cash_gap_expl(15))];
   disp(X);

    % compute gap in equity shares in model: baseline and all counterfactuals
    ucshare_gap = NaN(1,15);
    ucshare_gap(1)  = mean(ucsharecoh_sm)   -  mean(ucsharecoh_sf); 
    ucshare_gap(2)  = mean(ucshare_sm_dc1)  -  mean(ucshare_sf_dc1);
    ucshare_gap(3)  = mean(ucshare_sm_dc1b) -  mean(ucshare_sf_dc1b);
    ucshare_gap(4)  = mean(ucshare_sm_dc1c) -  mean(ucshare_sf_dc1c);
    ucshare_gap(5)  = mean(ucshare_sm_dc1d) -  mean(ucshare_sf_dc1d);
    ucshare_gap(6)  = mean(ucshare_sm_dc1e) -  mean(ucshare_sf_dc1e);
    ucshare_gap(7)  = mean(ucshare_sm_dc1f) -  mean(ucshare_sf_dc1f);
    ucshare_gap(8)  = mean(ucshare_sm_dc2)  -  mean(ucshare_sf_dc2);
    ucshare_gap(9)  = mean(ucshare_sm_dc3)  -  mean(ucshare_sf_dc3);
    ucshare_gap(10) = mean(ucshare_sm_dc4)  -  mean(ucshare_sf_dc4);
    ucshare_gap(11) = mean(ucshare_sm_dc5)  -  mean(ucshare_sf_dc5);
    ucshare_gap(12) = mean(ucshare_sm_dc6)  -  mean(ucshare_sf_dc6);
    ucshare_gap(13) = mean(ucshare_sm_dc7)  -  mean(ucshare_sf_dc7);
    ucshare_gap(14) = mean(ucshare_sm_dc8)  -  mean(ucshare_sf_dc8);
    ucshare_gap(15) = mean(ucshare_sm_dc9)  -  mean(ucshare_sf_dc9);
  
    X= ['Decomp share:   ','baseline: ' , num2str(ucshare_gap(1)),' income level: ' , ...
           num2str(ucshare_gap(2)), '  fixed term: ' , num2str(ucshare_gap(3)), ' age term: ' , num2str(ucshare_gap(4)), ...
         ' educ term: ' , num2str(ucshare_gap(5)),' only level: ' , num2str(ucshare_gap(6)), ...
         ' only slope: ' , num2str(ucshare_gap(7)),' income risk: ' , num2str(ucshare_gap(8)),...
         ' HH size: ' , num2str(ucshare_gap(9)),' marriage prob: ' , num2str(ucshare_gap(10)),... 
         ' marriage market: ' , num2str(ucshare_gap(11)), ' initial wealth: ' , num2str(ucshare_gap(12)),...
         ' education distribution: ' , num2str(ucshare_gap(13)),' medical expenses: ' , num2str(ucshare_gap(14)),...
         ' survival risk: ' , num2str(ucshare_gap(15))];
    disp(X);

    % what part of gap can be explained by individual channels?
    ucshare_gap_expl = NaN(1,15);
    ucshare_gap_expl(1)  = 1-(ucshare_gap(1)/ucshare_gap(1));
    ucshare_gap_expl(2)  = 1-(ucshare_gap(2)/ucshare_gap(1));
    ucshare_gap_expl(3)  = 1-(ucshare_gap(3)/ucshare_gap(1));
    ucshare_gap_expl(4)  = 1-(ucshare_gap(4)/ucshare_gap(1));
    ucshare_gap_expl(5)  = 1-(ucshare_gap(5)/ucshare_gap(1));
    ucshare_gap_expl(6)  = 1-(ucshare_gap(6)/ucshare_gap(1));
    ucshare_gap_expl(7)  = 1-(ucshare_gap(7)/ucshare_gap(1));
    ucshare_gap_expl(8)  = 1-(ucshare_gap(8)/ucshare_gap(1));
    ucshare_gap_expl(9)  = 1-(ucshare_gap(9)/ucshare_gap(1));
    ucshare_gap_expl(10) = 1-(ucshare_gap(10)/ucshare_gap(1));
    ucshare_gap_expl(11) = 1-(ucshare_gap(11)/ucshare_gap(1));
    ucshare_gap_expl(12) = 1-(ucshare_gap(12)/ucshare_gap(1));
    ucshare_gap_expl(13) = 1-(ucshare_gap(13)/ucshare_gap(1));
    ucshare_gap_expl(14) = 1-(ucshare_gap(14)/ucshare_gap(1));
    ucshare_gap_expl(15) = 1-(ucshare_gap(15)/ucshare_gap(1));
    
    X= ['Decomp share explained:   ','baseline: ' , num2str(ucshare_gap_expl(1)),' income level: ' , ...
           num2str(ucshare_gap_expl(2)), '  fixed term: ' , num2str(ucshare_gap_expl(3)), ' age term: ' , num2str(ucshare_gap_expl(4)), ...
         ' educ term: ' , num2str(ucshare_gap_expl(5)),' only level: ' , num2str(ucshare_gap_expl(6)), ...
         ' only slope: ' , num2str(ucshare_gap_expl(7)),' income risk: ' , num2str(ucshare_gap_expl(8)),...
         ' HH size: ' , num2str(ucshare_gap_expl(9)),' marriage prob: ' , num2str(ucshare_gap_expl(10)),... 
         ' marriage market: ' , num2str(ucshare_gap_expl(11)), ' initial wealth: ' , num2str(ucshare_gap_expl(12)),...
         ' education distribution: ' , num2str(ucshare_gap_expl(13)),' medical expenses: ' , num2str(ucshare_gap_expl(14)),...
         ' survival risk: ' , num2str(ucshare_gap_expl(15))];
    disp(X);

    %% %%%%%%%%% Composition vs. Policy Effect (Section 6.2)  %%%%%%%%%%%%
         
%% Policy vs. Composition Effect: Aggregate (Figure 5)

    % STOCK MARKET PARTICIPATION RATE
    fid = fopen('polcomp_agg\cohort\SMPcoh_sf.txt');
    SMP_sf_pc  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('polcomp_agg\cohort\SMPcoh_sm.txt');
    SMP_sm_pc  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('polcomp_agg\cohort\SMPcoh_c.txt');
    SMP_c_pc   = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    % CONDITIONAL EQUITY SHARE
    fid = fopen('polcomp_agg\cohort\sharecoh_sf.txt');
    share_sf_pc  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('polcomp_agg\cohort\sharecoh_sm.txt');
    share_sm_pc  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('polcomp_agg\cohort\sharecoh_c.txt');
    share_c_pc  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    % UNCONDITIONAL EQUITY SHARE
    fid = fopen('polcomp_agg\cohort\ucsharecoh_sf.txt');
    ucshare_sf_pc  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('polcomp_agg\cohort\ucsharecoh_sm.txt');
    ucshare_sm_pc  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('polcomp_agg\cohort\ucsharecoh_c.txt');
    ucshare_c_pc  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    %%%% PLOT RESULTS %%%%

        % Figure 5a
        figure(51)
        hold on
        plot(age, 100*(ucsharecoh_sm - ucsharecoh_sf), 'Linewidth', 2.5,'Color', 'black' )
        plot(age, 100*(ucshare_sf_pc - ucsharecoh_sf), 'Linewidth', 2.5,'Color', 'black','LineStyle', ':')
        plot(age, 100*(ucsharecoh_sm - ucshare_sf_pc), 'Linewidth', 2.5,'Color', 'black','LineStyle', '--')
        hold off
        xlim([30 64])
        ylim([-10 20])
        xlabel('Age','Interpreter','LaTex','Fontsize',24)
        ylabel('\%-points','Interpreter','LaTex','Fontsize',24)
        set(gca,'XGrid','off','YGrid','on','Fontsize',20)
        saveas(gcf,[figurefolder '\ucshare_polcomp.png'])
        hold off

        % Figure 5b
        figure(52)
        hold on
        plot(age, 100*(SMPcoh_sm - SMPcoh_sf), 'Linewidth', 2.5,'Color', 'black' )
        plot(age, 100*(SMP_sf_pc - SMPcoh_sf), 'Linewidth', 2.5,'Color', 'black','LineStyle', ':')
        plot(age, 100*(SMPcoh_sm - SMP_sf_pc), 'Linewidth', 2.5,'Color', 'black','LineStyle', '--')
        hold off
        xlim([30 64])
        ylim([-20 40])
        xlabel('Age','Interpreter','LaTex','Fontsize',24)
        ylabel('\%-points','Interpreter','LaTex','Fontsize',24)
        set(gca,'XGrid','off','YGrid','on','Fontsize',20)
        saveas(gcf,[figurefolder '\SMP_polcomp.png'])
        hold off

        % Figure 5c
        figure(53)
        hold on
        plot(age, 100*(sharecoh_sm - sharecoh_sf), 'Linewidth', 2.5,'Color', 'black' )
        plot(age, 100*(share_sf_pc - sharecoh_sf), 'Linewidth', 2.5,'Color', 'black','LineStyle', ':')
        plot(age, 100*(sharecoh_sm - share_sf_pc), 'Linewidth', 2.5,'Color', 'black','LineStyle', '--')
        hold off
        xlim([30 64])
        ylim([-10 20])
        xlabel('Age','Interpreter','LaTex','Fontsize',24)
        ylabel('\%-points','Interpreter','LaTex','Fontsize',24)
        set(gca,'XGrid','off','YGrid','on','Fontsize',20)
        leg = legend('total effect','policy effect','composition effect','Interpreter','LaTex');
        set(leg,'Interpreter','LaTex','Fontsize',19,'Location','NorthEast','NumColumns',1);
        saveas(gcf,[figurefolder '\share_polcomp.png'])
        hold off

      % policy effect, relative importance until age 40
        mean((ucshare_sf_pc(1:12)-ucsharecoh_sf(1:12))./(ucsharecoh_sm(1:12) - ucsharecoh_sf(1:12)))

%% Policy effect -- det. income channel (Figures 6, 15 and 16)

    % SMP RATE
    fid = fopen('polcomp_detinc\cohort\SMPcoh_sf.txt');
    SMP_sf_pc_income  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('polcomp_detinc\cohort\SMPcoh_sm.txt');
    SMP_sm_pc_income  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('polcomp_detinc\cohort\SMPcoh_c.txt');
    SMP_c_pc_income   = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    % CONDITIONAL EQUITY SHARE
    fid = fopen('polcomp_detinc\cohort\sharecoh_sf.txt');
    share_sf_pc_income  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('polcomp_detinc\cohort\sharecoh_sm.txt');
    share_sm_pc_income  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('polcomp_detinc\cohort\sharecoh_c.txt');
    share_c_pc_income  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    % UNCONDITIONAL EQUITY SHARE
    fid = fopen('polcomp_detinc\cohort\ucsharecoh_sf.txt');
    ucshare_sf_pc_income  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('polcomp_detinc\cohort\ucsharecoh_sm.txt');
    ucshare_sm_pc_income  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('polcomp_detinc\cohort\ucsharecoh_c.txt');
    ucshare_c_pc_income  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    %%%% PLOT RESULTS %%%%
        
        % Figure 6a
        figure(61)
        hold on
        plot(age, 100*(ucsharecoh_sm - ucsharecoh_sf), 'Linewidth', 2.5,'Color', 'black' )
        plot(age, 100*(ucsharecoh_sm - ucshare_sf_pc), 'Linewidth', 2.5,'Color', 'black','LineStyle', '--')
        plot(age, 100*(ucshare_sf_pc - ucshare_sf_pc_income), 'Linewidth', 2.5,'Color', 'black','LineStyle', ':')
        plot(age, 100*(ucshare_sf_pc_income - ucsharecoh_sf), 'Linewidth', 2.5,'Color', 'black','LineStyle', '-.')
        hold off
        xlim([30 64])
        ylim([-10 30])
        xlabel('Age','Interpreter','LaTex','Fontsize',24)
        ylabel('\%-points','Interpreter','LaTex','Fontsize',24)
        set(gca,'XGrid','off','YGrid','on','Fontsize',20)
        leg = legend('total effect','composition','policy: other factors','policy: det. income', 'Interpreter','LaTex');
        set(leg,'Interpreter','LaTex','Fontsize',19,'Location','NorthEast','NumColumns',1);
        saveas(gcf,[figurefolder '\ucshare_polcomp_income.png'])
        hold off

        % Figure 15a
        figure(151)
        hold on
        plot(age, 100*(SMPcoh_sm - SMPcoh_sf), 'Linewidth', 2.5,'Color', 'black' )
        plot(age, 100*(SMP_sf_pc - SMP_sf_pc_income), 'Linewidth', 2.5,'Color', 'black','LineStyle', ':')
        plot(age, 100*(SMP_sf_pc_income - SMPcoh_sf), 'Linewidth', 2.5,'Color', 'black','LineStyle', '-.')
        hold off
        xlim([30 64])
        ylim([-20 40])
        xlabel('Age','Interpreter','LaTex','Fontsize',24)
        ylabel('\%-points','Interpreter','LaTex','Fontsize',24)
        set(gca,'XGrid','off','YGrid','on','Fontsize',20)
        leg = legend('total effect','policy: other factors','policy: det. income', 'Interpreter','LaTex');
        set(leg,'Interpreter','LaTex','Fontsize',19,'Location','NorthEast','NumColumns',1);
        saveas(gcf,[figurefolder '\SMP_polcomp_income.png'])
        hold off

        % Figure 16a
        figure(161)
        hold on
        plot(age, 100*(sharecoh_sm - sharecoh_sf), 'Linewidth', 2.5,'Color', 'black' )
        plot(age, 100*(share_sf_pc - share_sf_pc_income), 'Linewidth', 2.5,'Color', 'black','LineStyle', ':')
        plot(age, 100*(share_sf_pc_income - sharecoh_sf), 'Linewidth', 2.5,'Color', 'black','LineStyle', '-.')
        hold off
        xlim([30 64])
        ylim([-10 20])
        xlabel('Age','Interpreter','LaTex','Fontsize',24)
        ylabel('\%-points','Interpreter','LaTex','Fontsize',24)
        set(gca,'XGrid','off','YGrid','on','Fontsize',20)
        leg = legend('total effect','policy: other factors','policy: det. income', 'Interpreter','LaTex');
        set(leg,'Interpreter','LaTex','Fontsize',19,'Location','NorthEast','NumColumns',1);
        saveas(gcf,[figurefolder '\share_polcomp_income.png'])
        hold off


%% Policy effect -- HH size channel (Figures 6, 15 and 16)    

    % SMP RATE
    fid = fopen('polcomp_hhsize\cohort\SMPcoh_sf.txt');
    SMP_sf_pc_HH  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('polcomp_hhsize\cohort\SMPcoh_sm.txt');
    SMP_sm_pc_HH  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('polcomp_hhsize\cohort\SMPcoh_c.txt');
    SMP_c_pc_HH   = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    % CONDITIONAL EQUITY SHARE
    fid = fopen('polcomp_hhsize\cohort\sharecoh_sf.txt');
    share_sf_pc_HH  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('polcomp_hhsize\cohort\sharecoh_sm.txt');
    share_sm_pc_HH  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('polcomp_hhsize\cohort\sharecoh_c.txt');
    share_c_pc_HH  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    % UNCONDITIONAL EQUITY SHARE
    fid = fopen('polcomp_hhsize\cohort\ucsharecoh_sf.txt');
    ucshare_sf_pc_HH  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('polcomp_hhsize\cohort\ucsharecoh_sm.txt');
    ucshare_sm_pc_HH  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('polcomp_hhsize\cohort\ucsharecoh_c.txt');
    ucshare_c_pc_HH  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

 
    %%%% PLOT RESULTS %%%%

        % Figure 6b 
        figure(62)
        hold on
        plot(age, 100*(ucsharecoh_sm - ucsharecoh_sf), 'Linewidth', 2.5,'Color', 'black' )
        plot(age, 100*(ucsharecoh_sm - ucshare_sf_pc), 'Linewidth', 2.5,'Color', 'black','LineStyle', '--')
        plot(age, 100*(ucshare_sf_pc - ucshare_sf_pc_HH), 'Linewidth', 2.5,'Color', 'black','LineStyle', ':')
        plot(age, 100*(ucshare_sf_pc_HH- ucsharecoh_sf), 'Linewidth', 2.5,'Color', 'black','LineStyle', '-.')
        hold off
        xlim([30 64])
        ylim([-10 30])
        xlabel('Age','Interpreter','LaTex','Fontsize',24)
        ylabel('\%-points','Interpreter','LaTex','Fontsize',24)
        set(gca,'XGrid','off','YGrid','on','Fontsize',20)
        leg = legend('total effect','composition','policy: other factors','policy: HH sizes', 'Interpreter','LaTex');
        set(leg,'Interpreter','LaTex','Fontsize',19,'Location','NorthEast','NumColumns',1);
        saveas(gcf,[figurefolder '\ucshare_polcomp_hhsize.png'])
        hold off

        % Figure 15b 
        figure(152)
        hold on
        plot(age, 100*(SMPcoh_sm - SMPcoh_sf), 'Linewidth', 2.5,'Color', 'black' )
        plot(age, 100*(SMP_sf_pc - SMP_sf_pc_HH), 'Linewidth', 2.5,'Color', 'black','LineStyle', ':')
        plot(age, 100*(SMP_sf_pc_HH- SMPcoh_sf), 'Linewidth', 2.5,'Color', 'black','LineStyle', '-.')
        hold off
        xlim([30 64])
        ylim([-20 40])
        xlabel('Age','Interpreter','LaTex','Fontsize',24)
        ylabel('\%-points','Interpreter','LaTex','Fontsize',24)
        set(gca,'XGrid','off','YGrid','on','Fontsize',20)
        leg = legend('total effect','policy: other factors','policy: HH sizes', 'Interpreter','LaTex');
        set(leg,'Interpreter','LaTex','Fontsize',19,'Location','NorthEast','NumColumns',1);
        saveas(gcf,[figurefolder '\SMP_polcomp_hhsize.png'])
        hold off

        % Figure 15c 
        figure(153)
        hold on
        plot(age, 100*(sharecoh_sm - sharecoh_sf), 'Linewidth', 2.5,'Color', 'black' )
        plot(age, 100*(share_sf_pc - share_sf_pc_HH), 'Linewidth', 2.5,'Color', 'black','LineStyle', ':')
        plot(age, 100*(share_sf_pc_HH - sharecoh_sf), 'Linewidth', 2.5,'Color', 'black','LineStyle', '-.')
        hold off
        xlim([30 64])
        ylim([-10 20])
        xlabel('Age','Interpreter','LaTex','Fontsize',24)
        ylabel('\%-points','Interpreter','LaTex','Fontsize',24)
        set(gca,'XGrid','off','YGrid','on','Fontsize',20)
        leg = legend('total effect','policy: other factors','policy: HH sizes', 'Interpreter','LaTex');
        set(leg,'Interpreter','LaTex','Fontsize',19,'Location','NorthEast','NumColumns',1);
        saveas(gcf,[figurefolder '\share_polcomp_hhsize.png'])
        hold off

%% Policy effect -- stoch. income channel (Figures 6, 15 and 16) 

    % SMP RATE
    fid = fopen('polcomp_risk\cohort\SMPcoh_sf.txt');
    SMP_sf_pc_risk  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('polcomp_risk\cohort\SMPcoh_sm.txt');
    SMP_sm_pc_risk  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('polcomp_risk\cohort\SMPcoh_c.txt');
    SMP_c_pc_risk   = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    % CONDITIONAL EQUITY SHARE
    fid = fopen('polcomp_risk\cohort\sharecoh_sf.txt');
    share_sf_pc_risk  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('polcomp_risk\cohort\sharecoh_sm.txt');
    share_sm_pc_risk  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('polcomp_risk\cohort\sharecoh_c.txt');
    share_c_pc_risk  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    % UNCONDITIONAL EQUITY SHARE
    fid = fopen('polcomp_risk\cohort\ucsharecoh_sf.txt');
    ucshare_sf_pc_risk  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('polcomp_risk\cohort\ucsharecoh_sm.txt');
    ucshare_sm_pc_risk  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    fid = fopen('polcomp_risk\cohort\ucsharecoh_c.txt');
    ucshare_c_pc_risk  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    %%%% PLOT RESULTS %%%%

        % Figure 6c
        figure(63)
        hold on
        plot(age, 100*(ucsharecoh_sm - ucsharecoh_sf), 'Linewidth', 2.5,'Color', 'black' )
        plot(age, 100*(ucsharecoh_sm - ucshare_sf_pc), 'Linewidth', 2.5,'Color', 'black','LineStyle', '--')
        plot(age, 100*(ucshare_sf_pc - ucshare_sf_pc_risk), 'Linewidth', 2.5,'Color', 'black','LineStyle', ':')
        plot(age, 100*(ucshare_sf_pc_risk - ucsharecoh_sf), 'Linewidth', 2.5,'Color', 'black','LineStyle', '-.')
        hold off
        xlim([30 64])
        ylim([-10 30])
        xlabel('Age','Interpreter','LaTex','Fontsize',24)
        ylabel('\%-points','Interpreter','LaTex','Fontsize',24)
        set(gca,'XGrid','off','YGrid','on','Fontsize',20)
        leg = legend('total effect','composition','policy: other factors','policy: stoch. income', 'Interpreter','LaTex');
        set(leg,'Interpreter','LaTex','Fontsize',19,'Location','NorthEast','NumColumns',1);
        saveas(gcf,[figurefolder '\ucshare_polcomp_risk.png'])
        hold off

        % Figure 15c
        figure(153)
        hold on
        plot(age, 100*(SMPcoh_sm - SMPcoh_sf), 'Linewidth', 2.5,'Color', 'black' )
        plot(age, 100*(SMP_sf_pc - SMP_sf_pc_risk), 'Linewidth', 2.5,'Color', 'black','LineStyle', ':')
        plot(age, 100*(SMP_sf_pc_risk - SMPcoh_sf), 'Linewidth', 2.5,'Color', 'black','LineStyle', '-.')
        hold off
        xlim([30 64])
        ylim([-20 40])
        xlabel('Age','Interpreter','LaTex','Fontsize',24)
        ylabel('\%-points','Interpreter','LaTex','Fontsize',24)
        set(gca,'XGrid','off','YGrid','on','Fontsize',20)
        leg = legend('total effect','policy: other factors','policy: stoch. income', 'Interpreter','LaTex');
        set(leg,'Interpreter','LaTex','Fontsize',19,'Location','NorthEast','NumColumns',1);
        saveas(gcf,[figurefolder '\SMP_polcomp_risk.png'])
        hold off

        % Figure 16c
        figure(163)
        hold on
        plot(age, 100*(sharecoh_sm - sharecoh_sf), 'Linewidth', 2.5,'Color', 'black' )
        plot(age, 100*(share_sf_pc - share_sf_pc_risk), 'Linewidth', 2.5,'Color', 'black','LineStyle', ':')
        plot(age, 100*(share_sf_pc_risk - sharecoh_sf), 'Linewidth', 2.5,'Color', 'black','LineStyle', '-.')
        hold off
        xlim([30 64])
        ylim([-10 20])
        xlabel('Age','Interpreter','LaTex','Fontsize',24)
        ylabel('\%-points','Interpreter','LaTex','Fontsize',24)
        set(gca,'XGrid','off','YGrid','on','Fontsize',20)
        leg = legend('total effect','policy: other factors','policy: stoch. income', 'Interpreter','LaTex');
        set(leg,'Interpreter','LaTex','Fontsize',19,'Location','NorthEast','NumColumns',1);
        saveas(gcf,[figurefolder '\share_polcomp_risk.png'])
        hold off
   
%% Simulate the dataset (to run regressions on simulated data - Table 4)   

 %%%% Export simulated data to .csv format to then run regressions in Stata %%%%

  % equity share
  fid = fopen('base\simdata\share_path_sf.txt', 'r');    
       share_path_sf  = cell2mat(textscan(fid, '%f')).';
       fclose(fid);  
       share_path_sf = reshape(share_path_sf,SS,JJ+1);
  csvwrite('share_path_sf.csv',share_path_sf,1,1);
 
  fid = fopen('base\simdata\share_path_sm.txt', 'r');    
       share_path_sm  = cell2mat(textscan(fid, '%f')).';
       fclose(fid);  
       share_path_sm = reshape(share_path_sm,SS,JJ+1);
  csvwrite('share_path_sm.csv',share_path_sm,1,1);
 
  fid = fopen('base\simdata\share_path_c.txt', 'r');    
       share_path_c  = cell2mat(textscan(fid, '%f')).';
       fclose(fid);  
       share_path_c = reshape(share_path_c,2*SS,JJ+1);
  csvwrite('share_path_c.csv',share_path_c,1,1);
    
  % assetgrid 
  fid = fopen('base\simdata\assetgrid.txt', 'r');    
       assetgrid = cell2mat(textscan(fid, '%f')).';
       fclose(fid);  
  csvwrite('assetgrid.csv',assetgrid,1,1);    
 
 
  % asset levels
  fid = fopen('base\simdata\grid_path_m.txt', 'r');    
        grid_path_m = cell2mat(textscan(fid, '%f')).';
        fclose(fid);  
  grid_path_m = reshape(grid_path_m,SS,JJ+1); 
  grid_path_m = grid_path_m*5000;
  csvwrite('grid_path_m.csv',grid_path_m,1,1); 
  
  fid = fopen('base\simdata\grid_path_f.txt', 'r');    
       grid_path_f = cell2mat(textscan(fid, '%f')).';
       fclose(fid);  
  grid_path_f = reshape(grid_path_f,SS,JJ+1);
  grid_path_f = grid_path_f*5000;
  csvwrite('grid_path_f.csv',grid_path_f,1,1); 
  
  fid = fopen('base\simdata\grid_path_c.txt', 'r');    
       grid_path_c = cell2mat(textscan(fid, '%f')).';
       fclose(fid); 
  grid_path_c = reshape(grid_path_c,2*SS,JJ+1); 
  grid_path_c = grid_path_c*5000;
  csvwrite('grid_path_c.csv',grid_path_c,1,1); 
   
  % stochastic income grid
  fid = fopen('base\afs_boom.txt', 'r');    
       afs_boom = cell2mat(textscan(fid, '%f')).';
       fclose(fid);  

  fid = fopen('base\afs_rec.txt', 'r');    
       afs_rec = cell2mat(textscan(fid, '%f')).';
       fclose(fid); 
 
  fid = fopen('base\ams_boom.txt', 'r');    
       ams_boom = cell2mat(textscan(fid, '%f')).';
       fclose(fid);  

  fid = fopen('base\ams_rec.txt', 'r');    
       ams_rec = cell2mat(textscan(fid, '%f')).';
       fclose(fid); 
  
  % chains for stochastic income realization 
  fid = fopen('base\chain_prodf.txt', 'r');    
        chain_prodf = cell2mat(textscan(fid, '%f')).';
        fclose(fid);  
        chain_prodf = reshape(chain_prodf,SS,JJ+1); 

  fid = fopen('base\chain_prodm.txt', 'r');    
        chain_prodm = cell2mat(textscan(fid, '%f')).';
        fclose(fid);  
        chain_prodm = reshape(chain_prodm,SS,JJ+1);

  % aggregrate state
  fid = fopen('base\chain_agg.txt', 'r');
        chain_agg  = cell2mat(textscan(fid, '%f')).';
        fclose(fid);  
        chain_agg = reshape(chain_agg,SS,JJ+1);

  % transform stoch. income indices to levels, dep. on aggregate state
    chain_sincf = zeros(SS,JJ);
    for t = 1: JJ+1
    for s = 1:SS
    if     (chain_agg(s,t) == 1) 
      chain_sincf(s,t) = afs_boom(chain_prodf(s,t));
    elseif (chain_agg(s,t) == 2) 
      chain_sincf(s,t) = afs_rec(chain_prodf(s,t));
    end 
    end  
    end  
    csvwrite('chain_prodf.csv',chain_sincf,1,1);
    
    chain_sincm = zeros(SS,JJ);
    for t = 1: JJ+1
    for s = 1:SS
    if     (chain_agg(s,t) == 1) 
       chain_sincm(s,t) = ams_boom(chain_prodm(s,t));
    elseif (chain_agg(s,t) == 2) 
       chain_sincm(s,t) = ams_rec(chain_prodm(s,t));
    end 
    end  
    end  
    csvwrite('chain_prodm.csv',chain_sincm,1,1);

%% Robustness Checks (Section 7) 

 % ROBUSTNESS 1: NO ASSET DROP IF ONE SPOUSE DIES

    % UNCONDITIONAL EQUITY SHARE
    fid = fopen('robustness_noassetdrop\cohort\ucsharecoh_sf.txt');
    ucshare_sf_rob1  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('robustness_noassetdrop\cohort\ucsharecoh_sm.txt');
    ucshare_sm_rob1  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    % ASSET HOLDINGS
    fid = fopen('robustness_noassetdrop\cohort\cashcoh_sf.txt');
    cash_sf_rob1  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
        
    fid = fopen('robustness_noassetdrop\cohort\cashcoh_sm.txt');
    cash_sm_rob1  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);


  % ROBUSTNESS 2: PARAMETERS OF THE BEQUEST MOTIVE
    % UNCONDITIONAL EQUITY SHARE
    fid = fopen('robustness_bequestparam\cohort\ucsharecoh_sf.txt');
    ucshare_sf_rob2  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('robustness_bequestparam\cohort\ucsharecoh_sm.txt');
    ucshare_sm_rob2  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    % ASSET HOLDINGS
    fid = fopen('robustness_bequestparam\cohort\cashcoh_sf.txt');
    cash_sf_rob2  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
        
    fid = fopen('robustness_bequestparam\cohort\cashcoh_sm.txt');
    cash_sm_rob2  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);


 % ROBUSTNESS 3: CORRELATION OF INCOME PROCESS WITHIN COUPLES

    % UNCONDITIONAL EQUITY SHARE
    fid = fopen('robustness_correlation\cohort\ucsharecoh_sf.txt');
    ucshare_sf_rob3  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('robustness_correlation\cohort\ucsharecoh_sm.txt');
    ucshare_sm_rob3  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    % ASSET HOLDINGS
    fid = fopen('robustness_correlation\cohort\cashcoh_sf.txt');
    cash_sf_rob3  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
        
    fid = fopen('robustness_correlation\cohort\cashcoh_sm.txt');
    cash_sm_rob3  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

  % ROBUSTNESS 4: LUMPSUM PAYMENT AT AGE 55

    % UNCONDITIONAL EQUITY SHARE
    fid = fopen('robustness_lumpsum\cohort\ucsharecoh_sf.txt');
    ucshare_sf_rob4  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('robustness_lumpsum\cohort\ucsharecoh_sm.txt');
    ucshare_sm_rob4  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    % ASSET HOLDINGS
    fid = fopen('robustness_lumpsum\cohort\cashcoh_sf.txt');
    cash_sf_rob4  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
        
    fid = fopen('robustness_lumpsum\cohort\cashcoh_sm.txt');
    cash_sm_rob4  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);


  % ROBUSTNESS 5: MEDICAL EXPENSES OF COUPLES

    % UNCONDITIONAL EQUITY SHARE
    fid = fopen('robustness_medexpen\cohort\ucsharecoh_sf.txt');
    ucshare_sf_rob5  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('robustness_medexpen\cohort\ucsharecoh_sm.txt');
    ucshare_sm_rob5  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    % ASSET HOLDINGS
    fid = fopen('robustness_medexpen\cohort\cashcoh_sf.txt');
    cash_sf_rob5  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
        
    fid = fopen('robustness_medexpen\cohort\cashcoh_sm.txt');
    cash_sm_rob5  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

  % ROBUSTNESS 6: RETIREMENT INCOME OF COUPLES

    % UNCONDITIONAL EQUITY SHARE
    fid = fopen('robustness_pension\cohort\ucsharecoh_sf.txt');
    ucshare_sf_rob6  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('robustness_pension\cohort\ucsharecoh_sm.txt');
    ucshare_sm_rob6  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    % ASSET HOLDINGS
    fid = fopen('robustness_pension\cohort\cashcoh_sf.txt');
    cash_sf_rob6  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
        
    fid = fopen('robustness_pension\cohort\cashcoh_sm.txt');
    cash_sm_rob6  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

  % ROBUSTNESS 7: ASSET SPLITTING RULE UPON DIVORCE

    % UNCONDITIONAL EQUITY SHARE
    fid = fopen('robustness_splittingrule\cohort\ucsharecoh_sf.txt');
    ucshare_sf_rob7  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
    
    fid = fopen('robustness_splittingrule\cohort\ucsharecoh_sm.txt');
    ucshare_sm_rob7  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

    % ASSET HOLDINGS
    fid = fopen('robustness_splittingrule\cohort\cashcoh_sf.txt');
    cash_sf_rob7  = cell2mat(textscan(fid, '%f')).';
    fclose(fid); 
            
    fid = fopen('robustness_splittingrule\cohort\cashcoh_sm.txt');
    cash_sm_rob7  = cell2mat(textscan(fid, '%f')).';
    fclose(fid);

   %%%% DISPLAY RESULTS (Table 7) %%%%

    % gender gap in asset holdings
    cash_gap_robust = NaN(1,8);
    cash_gap_robust(1) = (mean(cashcoh_sm)   -  mean(cashcoh_sf))*5; 
    cash_gap_robust(2) = (mean(cash_sm_rob1) -  mean(cash_sf_rob1))*5;
    cash_gap_robust(3) = (mean(cash_sm_rob2) -  mean(cash_sf_rob2))*5;
    cash_gap_robust(4) = (mean(cash_sm_rob3) -  mean(cash_sf_rob3))*5;
    cash_gap_robust(5) = (mean(cash_sm_rob4) -  mean(cash_sf_rob4))*5;
    cash_gap_robust(6) = (mean(cash_sm_rob5) -  mean(cash_sf_rob5))*5;
    cash_gap_robust(7) = (mean(cash_sm_rob6) -  mean(cash_sf_rob6))*5;
    cash_gap_robust(8) = (mean(cash_sm_rob7) -  mean(cash_sf_rob7))*5;

    % gender gap in equity shares
    ucshare_gap_robust = NaN(1,8);
    ucshare_gap_robust(1) = (mean(ucsharecoh_sm)   -  mean(ucsharecoh_sf)); 
    ucshare_gap_robust(2) = (mean(ucshare_sm_rob1) -  mean(ucshare_sf_rob1));
    ucshare_gap_robust(3) = (mean(ucshare_sm_rob2) -  mean(ucshare_sf_rob2));
    ucshare_gap_robust(4) = (mean(ucshare_sm_rob3) -  mean(ucshare_sf_rob3));
    ucshare_gap_robust(5) = (mean(ucshare_sm_rob4) -  mean(ucshare_sf_rob4));
    ucshare_gap_robust(6) = (mean(ucshare_sm_rob5) -  mean(ucshare_sf_rob5));
    ucshare_gap_robust(7) = (mean(ucshare_sm_rob6) -  mean(ucshare_sf_rob6));
    ucshare_gap_robust(8) = (mean(ucshare_sm_rob7) -  mean(ucshare_sf_rob7));

